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We show that the time evolution operator of kicked quantum systems, although 
a full matrix of size N x N, can be diagonalized with the help of a new method 



I based on a suitable combination of Fast Fourier Transform and Lanczos algorithm 

(N 

in just 0(A^ InA^) operations. It allows the diagonalization of matrizes of sizes up 

CN '. 

I to ~ 10® going far beyond the possibilities of standard diagonalization techniques 

' which need O(N^) operations. We have applied this method to the kicked Harper 

(N ■ 

? ■ I. INTRODUCTION 

o . 

. Periodically kicked quantum systems have played a prominent role in studies of quantum 



model revealing its intricate spectral properties. 
PACS numbers: 02.60.Dc, 05.45. +b 



chaos for the analysis of signatures of classical chaos in quantum systems 0. Important 
examples are the kicked rotator, i.e., the quantum version of Chirikov's standard map 
the kicked top [0], and the kicked Harper model [Q. Their quantum time evolution operator 
allows a fast numerical iteration on wave packets just like in their classical analogs, which 
are maps that can be iterated very quickly. The study of kicked sytems has led to many 
discoveries, a very spectacular one being the suppression of classical chaotic diffusion due to 
quantum mechanical interference in the kicked rotator [|1|, which later on was understood 
by a mapping onto the Anderson model of disordered systems and onto the supersym- 
metric nonlinear a model for quasi one-dimensional wires [^]. Another example supporting 
the importance of kicked systems appeared in the study of quantum signatures of classical 
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chaos in quasiperiodic systems, in which the kicked Harper model 0,0-|T^ yielded surpris- 
ing phenomena [§,|12,|1^,|1^,|1^] , e.g., metal- insulator transitions tuned by classical chaos. In 
the future, kicked quantum systems most probably will be the first to be used for numer- 
ical studies of the quantum signatures of mixed systems having a hierarchical phase space 
structure at the boundary between regular and chaotic motion. For example, the semiclas- 



sically predicted and experimentally observed [^,0 fractal properties of conductance 
fluctuations wait to be verified by numerical quantum calculations. 



In this paper we will show how the Lanczos algorithm for diagonalizing matrices p2 



which is well suited for sparse matrices, can be used to exploit the advantages of kicked 
quantum systems even though their time evolution operator is a full matrix. The key 
point, which has been overlooked throughout twenty years of studying quantum chaos with 
kicked systems, is the following observation: The most time consuming step of the Lanczos 
algorithm, a matrix-vector multiplication, is analogous to a step of the time evolution for 
some wave packet. For kicked systems the time evolution is known to be efficiently performed 
by using the Fast Fourier Transform (FFT) algorithm |23|. Consequently, by combining the 
Lanczos algorithm and the FFT, diagonalization of kicked systems only takes of the order 
of A^^ In operations and for a matrix of size = 10^ needs about two days on a standard 
workstation. By using parallel computers matrices of size = 10^ should be diagonalizable. 

As a first application of this new method, we will study the quantum signatures of 
classical chaos in quasiperiodic systems. For the kicked Harper model (KHM) Oj^the critical 
line, Borgonovi and Shepelyansky [|l^] convincingly concluded from diagonalizing matrices 
of sizes up to A^ = 1775 that in the limit N ^ oo the spectrum, in addition to pure point 
and absolutely continuous components, can also have a singular continuous component. The 
Lanczos-FFT method allows us to easily diagonalize the KHM for the much larger size of 
A^ = 51536. We will show that the indications for a singular continuous component found in 
Ref. [|1^ do not persist for larger system sizes. We will thereby reveal the KHM's intricate 
spectral properties and demonstrate the applicability of our method. 
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II. THE LANCZOS-FFT METHOD 



The numerical advantage of the time evolution of a periodically kicked system, e.g., a 
one-dimensional Hamiltonian of the form 

H{t)=T{p) + V{x)J25{t-n), (1) 

n 

is due to the fact, that its time-evolution operator for one period factorizes into 

f/ = exp |-^T (p)} exp {-^V^ (a^)} • (2) 

Here, h is typically an effective Planck's constant, p is the momentum and x the position 
operator. The second factor describes the kick, whereas the first factor governs the free 
evolution between kicks. If U is periodic in x, then U can be represented as an infinite 
matrix depending on a Bloch phase 6^, the so-called quasi-momentum. If, in addition, U is 
periodic with period A^, i.e., Ui+n,i'+n = Ui^i', the time-evolution operator can be reduced to 
an X matrix depending on a second Bloch phase 6p. This is, e.g., the case for the kicked 
rotator and the kicked Harper model, when h takes the value 2'kM/N . In order to obtain 
the time evolution it will be applied iteratively to a wave packet described by a vector of 
length N . This application in general requires N'^ operations for a single time step, whereas 
the factorization of U [Eq. (H)] allows a much faster implementation: In either momentum 
or position representation, application of the first and second factor of f/, resp., amounts 
to a vector-vector multiplication using just operations. Going back and forth between 
position and momentum representation one uses the very efficient FFT algorithm, which 



needs of the order of A^ln A^ operations ||23|. A detailed description of this standard method 
is given in the appendix including its extension to cases where A^ is not a power of 2. Thus, 
for propagating a wave packet one period in time just O(A^lnA^) operations are needed, 
which allows the study of time evolutions in systems with sizes up to A^ ~ 10®. 

On the other hand, this dramatic advantage of kicked quantum systems when studying 
time evolutions has not yet been fully exploited for determining the eigenvalues of U . So 
far this advantage has been used in this context to perform the time evolution for zN steps 



taking just zN"^ hiN operations and to Fourier transform the resulting time series giving the 
unitary eigenvalues. Unfortunately, for large matrices this method fails for two reasons: (i) 
resolution: very often, the finite resolution 1 / {zN) will not resolve closeby eigenvalues. Even 
worse, for quasiperiodic quantum systems with fractal spectra of dimension D < 1, in which 
the typical level spacing scales as this requires at least O ^A^^+^Z-^ In A^^ operations, 

(ii) storage: if U has localized eigenf unctions, one needs to store the time evolution at 
0{N) sites and needs storage of 0{zN'^) bytes, which for = 10^ and a relatively small 
2; = 10 amounts to 10^ Gigabyte by far exceeding conventional storage capacities. Standard 
diagonalization methods need O (N^) operations and have so far restricted to values of 
only a few thousand. 

We will now describe explicitly the algorithm which makes use of the factorization of U 
[Eq. (0)] to obtain all eigenvalues in 0{N'^ \n.N) operations. To this end we start by recalling 



the standard Lanczos algorithm for a general N x N complex Hermitian matrix H [^]. It 
generates a sequence of tridiagonal real symmetric matrices 

ai (32 
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whose eigenvalues can be calculated very efficiently by standard routines and converge 
towards the eigenvalues of H. The calculation of a„ and /5„ is based on the following 
recursion relations: Let ipo,ipi G with ipQ = and ipi a normalized random vector. 
Then for n > 1 we have a„ = (tpnyH^Jn), Pn = {ipn^Hipn-i), and V'n+i = ^n+i/||V^n+i||, 
where V'n+i = Hipn — oinipn — Pn'ipn-i- If the calculation is carried out in exact arithmetic, 
then L]sf has the same eigenvalues as H. In fact, in this case L^r is related to H via a 
similarity transformation i.e., \E'~^if \& = L^r with \& = {ipi . . . iPn)- Due to the instabilty 
of the algorithm for finite precision arithmetic, n has to be larger than A^ in order to 
approximate all eigenvalues of H with sufficient accuracy. In this case, Ln will possess 
spurious eigenvalues not approximating any eigenvalue of H . These so-called ghosts can be 



detected efficiently using standard techniques |2^. Even though this algorithm is designed 



for Hermitian matrices it can be used to obtain the spectrum of a unitary matrix U by 
applying it to = i + ^'^'^ ^ ~ h ~ ^0 ' eigenvalue e*"^ of U 

the eigenvalues of and H~ are costu and sin a;, resp. 

The most time consuming step in each iteration of the Lanczos algorithm is the calcu- 
lation of Hipn which in general needs A^^ operations. If Hipn can be computed in fewer 
operations as is the case for sparse matrices, the algorithm will become faster than standard 
diagonalizaton techniques. As described above, for kicked systems in which f/ is a full ma- 
trix, one can make use of its factorization [Eq.(0)] and the FFT algorithm to calculate Uipn 
and thus H^ipn in just O {N In N) operations. Hence, by using the Lanczos-FFT method 
the eigenvalues of U can be calculated in O {N"^ In A^) operations, such that systems of sizes 
up to A^ ~ 10^ are treatable. 

Let us finally mention that the numerical workload may be further reduced if U obeys 
certain symmetries, i) If the N x N matrix U is symmetric, then = RelJ and H = ImU 
such that f/^ no longer occurs. This is, e.g., the case for the kicked Harper model if one 
uses the time evolution operator U for one period starting halfway between two kicks with 
6p = 0. The spectrum, of course, is independent of the chosen time when the one period 
time evolution operator starts, ii) If the eigenvalues e^'^ of U are symmetric with respect 
to cij = 0, it suffices to diagonalize H^. Each of these symmetries reduces the number of 
operations by a factor of 2. This completes the description of our method to determine the 
spectrum of a kicked quantum system in O {N"^ In A^) operations. 

III. THE SPECTRUM OF THE KICKED HARPER MODEL 



We will now apply the Lanczos-FFT method to the kicked Harper model (KHM) P,[7|-p!8 
By calculating the quasienergies of a system of size A^ ~ 50000 we will give strong evidence 
that a former result on the spectral properties of the KHM obtained by treating matrices 
of sizes up to A^ ~ 2000 cannot be maintained. Since kicked systems with a larger size can 
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up to now only be treated by the Lanczos-FFT method, this demonstrates the value of this 
algorithm. 

As a generalization of the Harper model [p5|-p7| the KHM is widely used to study the 
influence of a chaotic classical limit on the spectral properties of a quasiperiodic quantum 
system [^,|TT],|T2|,|I^|TB[ . The Harper model is described by a tight-binding Hamiltonian with 
potential Vn = V cos (2'Kan + u) at site n. The Hamiltonian of the KHM is given by 

H = L cosp + K cosx^^6 {t — n) , (4) 

n 

i.e., it is of the form (|l]) with T (p) = L cosp and V (x) = K cosx. For small values of K and 
L its spectrum is closely related to the spectrum of the Harper model 0. To be explicit, 
for small L/h and K/h the quasienergies of the KHM are given up to a factor L/2h by 
the eigenenergies of the Harper model with V = 2K/L, a = h/2TT and u = O^- Therefore, 
in the case of small K and L the spectral properties of the KHM are the same as for the 
Harper model: For a typical irrational /l/27r and fixed quasimomentum 9^^ the spectrum is 
pure point if > L, absolutely continuous ii L > and singular continuous for K = L 
with a fractal dimension given by that of the Harper model. Due to a duality property of 
the KHM pH] on the critical line K = L, the spectrum remains singular continuous for all 



values of K = L. Its fractal dimension, though, changes with increasing K ||TT|]. Off the 



critical line, for increasing values of K and L the spectral properties of the KHM differ from 
the Harper model in that in some regions pure point and absolutely continuous components 
coexist The changes in the spectral type for K ^ L can be understood with the 

help of avoided band crossings ( ABCs) |jl8[ . These are present in the spectrum of the KHM 



because its classical version generates chaotic dynamics (see Fig. 5 in Ref. 1|T8[). As has been 
shown in Ref. |T^ they may tune metal-insulator transitions, thus explaining the coexistence 
of pure point and absolutely continuous spectral components. Since the fractal dimensions 
of a spectrum may also change due to ABCs |T^, they provide a common explanation for 
all the changes in the quasienergy spectrum of the KHM mentioned above. 

Based on a numerical analysis of matrices of sizes up to = 1775, Borgonovi and 



Shepelyansky found indications that even a singular continuous component should exist 



off the critical line ||T^. In particular, they have studied the case of = 4, L = 7, 
h = 2n/(6 + (Tg), where ac = {V^ + l)/2 is the Golden Mean. For these parameters 
they have calculated the inverse participation ratio ^ = J2n \i^n\'^ of the eigenf unctions ip 
for the approximants 27r55/4 19 , 27rl44/1097 and 27r233/1775 of ^ = 27r/(6 + aa) yielding 
matrices of sizes N = 419, 1097 and 1775, resp. In some regions of the spectrum they 
have found the inverse participation ratio of the corresponding eigenf unctions to decrease 
as A^~^, indicating that there the spectrum is absolutely continuous, while in other regions 
the inverse participation ratio remains constant, indicating the spectrum to be pure point. 
Most interestingly, for quasienergies in the interval [2.4, 2.8] the inverse participation ratio 
of the corresponding eigenf unctions scaled with some exponent larger than -1 and smaller 
than 0, indicating a singular continuous spectral component in the limit — > oo. 

Using the Lanczos-FFT method we have calculated the quasienergy spectrum for the 
same values of K and L but for h = 27r6765/51536 that is a much higher approximant 
of 27r/(6 + (Jc)- This yields a matrix of size = 51536, by far exceeding the size of a 
matrix treatable by standard methods. In Fig. 1 we show the quasienergy spectrum versus 
the Bloch-phase Or^ and several magnifications of the quasienergy interval [2.4,2.8]. Note 
that the isolated dots present in the figure are ghosts produced by the Lanczos algorithm. 
From a plot of the eigenvalues of a periodic approximant as given here, one can infer the 
spectral type of the quasiperiodic system as follows: The spectrum is pure point if the width 
of Bloch bands is negligible as compared to their spacings. It is absolutely continuous if 
the spacings are negligible compared to the band widths. The self-similar structure of a 
singular continuous spectrum leads to band widths and spacings of comparable size as can 
be found, e.g., for the Harper model at the critical point. In Fig. 1 on intermediate scales the 
spectrum displays a self-similar structure in the quasienergy interval [2.4,2.8], confirming 
the findings of Ref. |]T^. As the final magnifications show, however, on smaller scales the self- 



similarity vanishes in many places and there the spectrum consists of either levels or bands 
corresponding to pure point and absolutely continuous components, resp. Note, however. 



that we also find small regions, where the self-similarity does not end even for our large 
matrix size (not shown in Fig. 1). This indicates that the fraction of self-similar scaling 
regions is decreasing with matrix size and questions the occurrence of a singular continuous 
component in the limit — > oo. How can one understand this surprising property of the 
KHM? 

This self-similar structure of the spectrum on intermediate scales can be understood with 
the help of avoided band crossings. As already mentioned, ABCs tune metal-insulator tran- 
sitions in the KHM leading to the coexistence of pure point and absolutely continuous com- 
ponents of the spectrum for 7^ L. Numerical investigations show that these components 
do not overlap. Hence, there will be a transistion manifold in the (quasienergy,ir,L)-space 
separating pure point from absolutely continuous regions. If this transition manifold, which 
may have a very complicated shape, would coincide with a line of fixed K and L over some 
quasienergy region, one would have a singular continuous component in this spectral region. 
Unless there is a special symmetry present (like for K — L), however, this is very unlikely 
to happen. Typically, a line of fixed K and L will have many intersections with this mani- 
fold, yielding quasienergy regions with either pure point or absolutely continuous spectrum. 
Thus there is no singular continuous component but the absolutely continuous and pure 
point components can be intermingled in a complicated manner depending on the shape of 
the transition manifold. The closer one is to an intersection with the transition manifold, the 
finer will be the scale where the self-similar looking spectrum decides for being pure point 
or absolutely continuous. For this reason the spectrum looks self-similar on intermediate 
scales, whereas on smaller scales pure point and absolutely continuous components show 
up. This is just what can be seen in the quasienergy spectrum shown in Fig. 1. This also 
explains why for some regions of the spectrum that are even closer to such an intersection, 
we find self-similar structure on all scales of our finite matrix. In the limit of an infinite 
matrix size the fraction of self-similar scaling regions should tend to zero and there should 
be no singular continuous component. 

In order to go beyond this qualitative understanding, we plan for the future to study 
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quantitatively how this fraction decays. Also, it will be interesting to map out the transition 
manifold in the (quasienergy,ir,L)-space by studying the spectrum for large matrices. One 
may wonder if it has a finite or an infinite number of intersections with a line of fixed K 
and L. 



IV. CONCLUDING REMARKS 

In this article we have presented a very powerful method for diagonalizing kicked quantum 
systems of sizes up to ~ 10^ on standard present day workstations. The necessity of 
diagonalizing matrices of this size was then demonstrated by the example of the KHM. Its 
intricate spectral structure was revealed and found to have its reason in the occurence of 
avoided band crossings. As was alluded to in the introduction, the possibility of obtaining 
the eigenvalues of large matrices is of great interest also in other areas of quantum chaos. 



Currently, our method is applied to the investigation of quantum fractal eigenstates 
In the future, some possible applications include the study of level statistics for very large 
spacings, the search for quantum signatures of the hierarchical phase space of mixed systems, 
and the study of predictions for disordered systems via the kicked rotator We express 
our hope that the Lanczos-FFT method will play a useful role in the field of quantum chaos. 

This research was carried out partially during two of the authors' stay at the Institute for 
Theoretical Physics, Santa Barbara. T. G. and R. K. gratefully acknowledge the hospitality 
of the ITP and its members. It was supported by the NSF under Grant No. PHY94-07194 
and in part by the Deutsche Forschungsgemeinschaft. 



V. APPENDIX 

In this appendix we give the explicit expressions usually employed for calculating Uip, 
when U may be factorized according to Eq. (j^). In, e.g., momentum representation we have 

t/ = exp |--^T Fp^.exp j-l^ (5) 

9 



with the effective Planck's constant H and Fp^x and F^^p denoting the Fourier transform 
from position to momentum space and vice versa, resp. If the first and the third factor 
on the r.h.s. are periodic with period A^, the operator reduces to a matrix of size N x N 
depending on two Bloch phases. Then, these factors are diagonal matrices with elements 
exp {-■ ir (^h{l + 1^))} and exp [-{V (27r{k + ^)/N^}, resp., in which k,l = 0, . . . , N -1, 
and Ox, Op G [0,27r] are the Bloch phases in position and momentum representation, resp. 
The Fourier transform Fp^x is the matrix 

= 1 (6) 



which factorizes in a scalar, in two diagonal matrices, and a standard Fourier transform, 
which can be implemented using the FFT algorithm. {Fx^p)^i is simply the conjugate 
complex of {Fp^x)ik such that in the final expression for U [Eq.(^] the first two terms of 
Eq.(0) cancel each other. 

The advantage of using the FFT algorithm exists only if is a power of 2; otherwise it 
is less efficient. In situations where N is not a power of 2, as is typically the case for the 
kicked Harper model where an irrational value of h/27T is approximated by a rational M/N, 
one may use the following strategy: One replaces the Fourier transforms in by Fourier 
transforms in where N > 2N is a power of 2. Before the first FFT one extends the 
vector of size to a vector of size by adding zeros and after the second FFT one retains 
only the first A^ elements of the resulting vector. In addition, one has to replace the diagonal 
matrix of size N x N with elements exp [-{V (2Ti{k + ^)/N^} by a diagonal matrix of 
size N X N with elements 



^ Eexp|-^V^f27rfm+^)/A^) le-^-^^'/V-^^'^/^, (8) 



in which k = 0, . . . ,N — 1. 
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FIGURES 




FIG. 1. Quasienergies u) vs. Bloch phase 9xN of the kicked Harper model [Eq.(^)] for = 4, 
L = 7, ?i = 27r 6765/51536, and ^^ = 0. Successive magnifications show fractal behaviour on several 
scales, whereas the last magnification reveals either localized eigenfunctions (no phase dependence 
of the quasienergies) or extended eigenfunctions (strong phase dependence and no visible gaps). 
The magnifications are for the intervals [2.56,2.76], [2.594,2.61], [2.6076,2.6092], [2.60874,2.60891], 
[2.59425, 2.59585] and [2.59555, 2.5958]. The isolated points are ghosts due to the instability of the 
Lanczos algorithm. 
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